#ifndef _EQUATION_SOLVER_H_
#define _EQUATION_SOLVER_H_

#include <iostream>
#include <cmath>

#define maximal 1e10
const static double epsl = 1e-8;

class Function{//虚拟函数类
public:
    virtual double operator () (const double &x) const = 0;
    virtual double posdev( const double &x) const{
        return ((*this)(x+epsl)-(*this)(x))/epsl;}
    virtual double negdev( const double &x) const{
        return ((*this)(x)-(*this)(x-epsl))/epsl;}
};

class zero: public Function{//零函数，用于测试
public:
    double operator () (const double &x) const {return 0;}
    double posdev( const double &x) const {return 0;}
    double negdev( const double &x) const {return 0;}
};

class EquationSolver{
protected:
    double e;//函数精度
    int iteration_Max;
    Function &f;
public:
    virtual double root() = 0;
    EquationSolver( Function &f, int Max = 1000, double e = epsl):
    f(f), e(e), iteration_Max(Max) {};
};

class BisectionSolver : public EquationSolver{
private:
    double a, b, d;//默认搜索区间[0,1] ，d为最小搜索区间宽度
public:
    BisectionSolver( Function &f, double a = 0, double b = 1, int Max = 1000, double e = epsl, double d = epsl):
    EquationSolver(f, Max, e), a(a), b(b), d(d) {};
    bool check(){
        if( f(a)*f(b) >= 0){
            std::cerr <<  "Invalid input! Programe will return a meaningless value." << std::endl;
            return false;
        }
        return true;
    }
    double root(){//如果区间内含有多个零点则返回最左侧的零点
        bool sign = check();
        int count = 0;
        double c, l = f(a), m, r = f(b);
        do{
            c = (a+b)/2;
            m = f(c);
            if( l * m < 0) b = c, r = m;
            else a = c, l = m; 
        }while( ++count <= iteration_Max && fabs(m) > e && fabs(b-a) > d);
        if(sign){
            std::cerr << "Result with bisection method:" << std::endl;
            std::cerr << "Iteration step:" <<  count << std::endl;
            std::cerr << "Root:" << c << std::endl;
            std::cerr << "Value:" << m << std::endl; 
        }
        return c;
    }
};

class NewtonSolver: protected EquationSolver{
private:
    double x_0, d;//默认初值为0
public:
    NewtonSolver( Function &f, double x = 0, int Max = 1000, double e = epsl, double d = epsl):
    EquationSolver(f, Max, e), x_0(x), d(d) {};
    double check( double &x){
        //std::cerr<<"check"<<std::endl;
        if( fabs(f.posdev(x)-f.negdev(x)) > 10000*e){
            std::cerr <<  "Non-derivable Function! Programe will return a meaningless value." << std::endl;
            return 0;
        }
        double d = (f.posdev(x)+f.negdev(x))/2;
        if( d == 0){
            std::cerr <<  "Not smooth Function! Programe will return a meaningless value." << std::endl;
        }
        //std::cerr<<"checked"<<std::endl;
        return d;
    }
    double deriv( double &x){
        //std::cerr<<"deriv"<<std::endl;
        if( (d=check(x)) != 0){
            return d;
        }
        return maximal;//返回无穷大
    }
    double root(){
        int count = 0;
        double x = x_0, value, dval;
        do{
            value = f(x);
            dval = deriv(x);
            x = x - value/dval;
        }while( ++count < iteration_Max && fabs(value) > e && fabs(value/dval) > d);
        if( dval < maximal){
            std::cerr << "Result with Newton method:" << std::endl;
            std::cerr << "Iteration step:" <<  count << std::endl;
            std::cerr << "Root:" << x << std::endl;
            std::cerr << "Value:" << value << std::endl;
        }
        return x;
    }
};

class SecantSolver: protected EquationSolver{
private:
    double a = 0, b = 1, d = epsl;
public:
    SecantSolver( Function &f, double a = 0, double b = 1,int Max = 1000, double e = epsl, double d = epsl):
    EquationSolver(f, Max, e), a(a), b(b), d(d) {};
    double root(){
        int count = 0;
        double high = f(a), low = f(b), dis;
        while( ++count < iteration_Max && fabs(low) > e && fabs(dis=(a-b)) > d){
            if( fabs(low) > fabs(high) ){
                std::swap(a, b);
                std::swap(low, high);
            }
            a = b;
            b = b - low*dis/(low-high);
            high = low;
            low = f(b);
        }
        std::cerr << "Result with Secant method:" << std::endl;
        std::cerr << "Iteration step:" <<  count << std::endl;
        std::cerr << "Root:" << b << std::endl;
        std::cerr << "Value:" << low << std::endl;
        return b;
    }
};


#endif